Ab-initio molecular dynamics for high-pressure liquid Hydrogen 
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We introduce an efficient scheme for the molecular dynamics of electronic systems by means of 
quantum Monte Carlo. The evaluation of the (Born-Oppenheimer) forces acting on the ionic posi- 
tions is achieved by two main ingredients: i) the forces are computed with finite and small variance, 
which allows the simulation of a a large number of atoms, ii) the statistical noise corresponding to 
the forces is used to drive the dynamics at finite temperature by means of an appropriate friction 
matrix. A first application to the high-density phase of Hydrogen is given, supporting the stability 
of the liquid phase at ~ SOOGPa and ~ 400K. 



O 



D 

i 

CO 



I 

C 

O 

o 



> 
o 
o 

oo 
m 
o 

o 

03 



i 

C 

O 

o 



X 



PACS numbers: 47.11.Mn, 02.70.Ss, 61.20.Ja, 62.50.+p 

The phase diagram of Hydrogen at high pressure is 
still under intense study from the experimental and the- 
oretical point of view. In particular in the low tempera- 
ture high-pressure regime there is yet no clear evidence 
of a metallic atomic solid, and a suggestion was given 
in Ref.(|2|) that the liquid phase is instead more stable. 
Indeed, for high pressures around 300GPa, a two fluid 
(proton and electron) superconducting phase of the con- 
ventional type, namely induced by the strong electron- 
phonon coupling, has been conjectured^]. In this work 
we use an improved ab-initio molecular dynamics (AMD) 
by using accurate forces computed by Quantum Monte 
Carlo (QMC). We present preliminary results, showing 
that the liquid phase is energetically more stable, due 
to the strong electron correlation, at least within the 
Resonating Valence Bond (RVB) variational approach[l|, 
which is very accurate also in the solid phase. 

AMD is well established as a powerful tool to investi- 
gate many-body condensed matter systems. Indeed, pre- 
vious attempts to apply Quantum Monte Carlo (QMC) 
for the dynamics of ions [J] or for their thermodynamic 
properties^] are known, but they were limited to small 
number N of electrons or to total energy corrections of 
the AMD trajectories, namely without the explicit cal- 
culation of the forces. 

Calculation of forces with finite variance. The simplest 
method for accurate calculations within QMC, is given by 
the so called variational Monte Carlo (VMC), which al- 
lows to compute the variational energy expectation value 



VMC 



(■iPt\4>t) 



of a highly accurate correlated wave 
function (WF) ipx by means of a statistical approach: 
electronic configurations {x}, with given electron posi- 
tions ft and spins o~i = ±1/2 for i = 1, • • • N, are usually 
generated by the Metropolis algorithm according to the 
probability density \i x oc iPt(x) 2 . Then Evmc is com- 
puted by averaging statistically over /i x the so called local 
energy e L (x) = , namely E VMC = J dfi x e L (x), 

where J d[i x indicates conventionally the 3iV multidimen- 
sional integral over the electronic coordinates weighted 



by i/jj,( x )- hi the present work we assume that the WF 
iPt{x) = (x\iPt) = J x det A is given by a correlated Jas- 
trow factor J times a determinant D of a N x N matrix 
A, such as for instance a Slater determinant. The main 
ideas of this approach can be straightforwardly general- 
ized to more complicated and more accurate WF's. 
The efficient calculation of the energy derivatives, 



namely the forces fs. 
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where Na is the number of atoms, is the most impor- 
tant ingredient for the AMD. Within VMC they can be 
computed by simple differentiation of Evmc, using that 
not only the Hamiltonian H but also ipT depend explic- 
itly on the atomic positions Ri. This leads to two dif- 



ferent contributions to the force = /Jf 



HF 



+ /£, the 



Hellmann-Feynman / 



HF 



and the Pulay one f? , where: 



fHF 



/J = -2 J d^ x (e L (x) - E V Mc)d Ai log\Mx)\(2) 

However in order to obtain a statistically meaningful 
average, namely with finite variance, some manipula- 
tions are necessary because the first integrand may di- 
verge when the atoms are close to some electronic posi- 
tions, whereas the second integrand is analogously un- 
bounded when a configuration x approaches the nodal 
surface determined by iPt( x ) — 0. By denning with d 
(8) the distance of x from the nodal region (the mini- 
mum electron-atom distance), eL(x),dslogipT{x) —1/d 
( (x\dfi.H\x) ~ l/<5 2 ), whereas fi x ~ d 2 (fi x ~ <5 2 ), lead- 
ing to an unbounded integral of the square integrand in 
Eq.© (EqfT]), namely to infinite variance. The infinite 
variance problem in Eq.([T]) was solved in several ways. 
Here we adopt a very elegant and efficient scheme pro- 
posed by Caffarel and AssarafQ. Instead the infinite 
variance problem in Eq.([2]) was not considered so far, 
and this is clearly a problem for a meaningful definition 
of ionic AMD consistent with QMC forces. 

In this letter we solve this problem in the following 
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simple way, by using the so called re-weighting method. 
We use a different probability distribution oc ipaix) 2 , 
determined by a guiding function ipa(x): 



iP G (x)=R*(x)(^ T (x)/R(x)) 



(3) 



where i?(a;) oc iPt(x) — > for d — > is a "measure" 
of the distance from the nodal surface ipr(x) = 0. By 
assumption ^ may vanish only when det A — ( J > 0) 
and therefore i?(x) is chosen to depend only on A. For 
reasons that will become clear later on we have adopted 
the following expression: 



R(x) = 1/ 



\ 
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(4) 



Then the guiding function is defined by properly regular- 
izing R(x), namely: 



R e (x) 



R(x) if R(x) > e 
e(R(x)/e) R W ( if R(x) < e 



(5) 



The non obvious regularization for R(x) < e instead of 
e.g. R e (x) = Max[e, R(x)] was considered in order to sat- 
isfy the continuity of the first derivative of i^c {%) when 
R(x) — e, thus ensuring that iPg(%) remains as close 
as possible to the trial function ipx- In this way the 
Metropolis algorithm can be applied for generating con- 
figurations according to a slightly different probability 
(i e (x) and the exact expression of f? can be obtained by 
the so called umbrella average: 



J Ri 



-2 / d/4 S(x)(e L (x) - E VM c)d^logip T {x 
M4 W) 



(6) 



Now, the re-weighting factor S(x) = (iPt(x)/4'g(x))' 2 = 
Min[l,(R(x)/e)] 2 ~ 2R ^ x ' > ' e oc d 2 , cancels out the diver- 
gence of the integrand, that was instead present in Eq.0. 
Hence the mentioned integrands in the numerator and 
S(x) ( < 1) in the denominator of Eq.© represent 
bounded random variables and have obviously finite vari- 
ance. In this way the problem of infinite variance is 
definitely solved within this simple re-weighting scheme. 
Moreover, in the present method R(x) is not related to 
an overall factor of the total WF, such as the total de- 
terminant, defined in a very wide range of values ~ e"" 
over the various configurations. It is instead obtained 
by using a quantity R(x) ~ -k with small fluctuations. 
Therefore the present scheme is particularly efficient and 
stable also for large N. 

We show in Fig.fl} a comparison of several methods 
for computing the Pulay force component acting on a 
Hydrogen proton at r s — 1.31 in a bec lattice. As it is 
clear in the plot for the N = 128 case, the difference be- 
tween a method with finite variance and the standard one 
with infinite variance is evident. Moreover for N = 250 



the simpler choice R(x) = \detA\ with finite (but large) 
variance is clearly very inefficient due to the difficulty to 
cross from the region with R(x) < e, where the integrand 
almost vanishes (see Fig[T]), to the one with R(x) > e and 
viceversa. Instead in the present scheme an appropriate 
choice of e, such that J d/i x S(x) ~ 1/2, allows frequent 
barrier crossings any few Metropolis steps also for large 
N. In this way the 3Na X 3Na correlation matrix ckqmc, 
defining the statistical correlation between the force com- 
ponents, can be efficiently evaluated: 

a QM c(R) =< (fg.- < fg. < fn, >) > (7) 

where the brackets <> indicate the statistical average 
over the QMC samples. The correlation matrix ckqmc, 
that within the conventional method is not even defined, 
will be a fundamental ingredient for a consistent AMD 
with QMC forces and therefore the solution of the infi- 
nite variance problem is particularly important for this 
purpose. 
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N=128 



a R(x)=l 

R(x) from Eq.5 




□ R(x)=| det A| 
R(i) from Eq.5 
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FIG. 1: (color online). Evolution of the integrand in Eq.© 
as a function of the Monte Carlo iterations. Each new sample 
is obtained after 2iV Metropolis trials. 



Langevin dynamics. In the following derivation we as- 
sume that ions have unit mass, that can be generally ob- 
tained by e.g. a simple rescaling of lengths for each ion 
independently. For clarity and compactness of notations, 
we also omit the ionic subindices i when not explicitly 
necessary. Moreover matrices (vectors) are indicated by 
a bar (arrow) over the corresponding symbols, and the 
matrix-vector product is also implicitly understood. We 
start therefore by the following AMD equations for the 
ion coordinates R and velocities v: 



= ~j(R)v + f(R)+Tj(t) 



R 



(8) 
(9) 



By using the fluctuation-dissipation theorem the friction 
matrix 7 is related to the temperature T (henceforth the 
Boltzmann constant fee = 1) by the relation: 

i(R) = ^"(i?) (10) 
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where a(R) is generally a symmetric correlation matrix: 

<ff i (t)ff j (t')>=S(t-t')a(R). (11) 

It is important to emphasize that, as a remarkable gener- 
alization of the standard AMD used in , m the present 
approach the friction matrix 7, may depend explicitly on 
the ion positions R, so that Eq. (fTU|) can be satisfied even 
for a generic correlation matrix a(R). In fact, since forces 
are computed by QMC, we can assume that there exists 
also a QMC contribution to a(R): 



a (R) = a + A olqmc (R) 



(12) 



where Ao > and ao is the identity matrix / up to 
another positive constant ao, ao = olqI. 

Integration of the Langevin dynamics. In the interval 
t n - A* < t < t n + A*, for At small, the positions R 
are changing a little and, within a good approximation, 
we can neglect the R dependence in the RHS of Eq.©. 
Moreover the velocities v n are computed at half-integer 
times t n — whereas coordinates R n are assumed to be 
defined at integer times R n — R(t n ). Then the integra- 
tion of Eq. ijHJ in the mentioned intervals can be recasted 
in the following useful form, where the force components 
are corrected by appropriate noisy vectors: 



V n +1 
Rn+1 
f 



(13) 
(14) 
(15) 

(16) 



By using that a = 2T7 from Eq. flOl) . the correlator defin- 
ing the discrete (time integrated) noise ff can be com- 
puted explicitly and is given by: 



= e^ At v n +f(f(R n )+^) 
= R n + Atv n+1 + 0(A 3 ) 
= ^(l-e-^*) 

2sinh(A* 7 ) 



< fjifjj >= 2T 7 



2 sinh(At 7 ) 
4sinh(A*7) 2 



(17) 



This means that the QMC noise has to be corrected in a 
non trivial way as explained in the following. 

Noise correction. The QMC noise is given during the 
simulation, and therefore in order to follow the correct 
dynamics another noise ff ext has to be added to the noisy 
force components in a way that the total integrated noise 
is the correct expression flTTl) . i.e. fj = ff ext + ffQuc- By 
using that the QMC noise in Eq.© is obviously inde- 
pendent of the external noise, we easily obtain the corre- 
sponding correlation matrix: 



< rfi Xt Vf xt >- 



OLQMC 



(18) 



On the other hand, after substituting the expression ([12]) 
in Eq. pTJ)) 7 = («o + Ao o:qmc) and using the ex- 
pression (|17p for a', we finally obtain a positive definite 



TABLE I: Comparison of the total energy per proton 
(Hartree) for Hydrogen in the bcc lattice at r s = 1.31 com- 
pared with the published ones with lowest energy (to our 
knowledge). All energies are in Hartree. 



N 


Evmc/Na 


Evmc/N a \12] 


Edmc/Na 


E DM c/Na\12] 


16 
54 
128 
250 


-0.48875(5) 
-0.53573(2) 
-0.49495(1) 
-0.49740(2) 


-0.4878(1) 
-0.5353(2) 
-0.4947(2) 


-0.49164(4) 
-0.53805(4) 
-0.49661(3) 
-0.49923(2) 


-0.4905(1) 
-0.5390(5) 
-0.4978(4) 



matrix in Eq. (TT8")) for At < Ao.fH Hence ff ext is a generic 
Gaussian correlated noise that can be easily sampled by 
standard algorithms. After that the random vector ff ext 
is added to the force / + f]QMC obtained by QMC, and 
replaces f + fj in Eq. (fl~3"|) . This finally allows to obtain 
an accurate AMD with a corresponding small time step 
error. In our approach the choice ao = is also allowed 
but, following Refs.(|7||9|), much better performances of 
the AMD are obtained with non zero ao > and/or 
Ao > At, namely with an external noise larger than the 
smallest possible one (ao = 0, Ao = At). The main 
advantage of this technique is that, at each iteration, by 
means of Eq. fTO]) . the statistical noise on the total energy 
and forces (see Figf2]) can be much larger than the tar- 
get temperature T, and this allows to improve the QMC 
efficiency by several orders of magnitude. 

Optimization of the WF. In the following examples we 
use a variational WF J x det A that is able to provide a 
very accurate description of the correlation energy, due to 
a particularly efficient choice of the determinant factor, 
that allows to describe the RVB correlations. fiol [ill] The 
WF contains several variational parameters, indicated by 
a vector /3, that have to be consistently optimized during 
the AMD. The Jastrow factor J used here depends both 
on the charge and spin densities, for a total of ~ N\ 
variational parameters. We employ periodic boundary 
conditions; for each proton we use two periodic gaussians 
centered at each ionic position, whereas for the Jastrow 
we use only one Gaussian. As it is shown in the table, the 
accuracy of our WF is remarkable. Indeed the small dif- 
ference between the so called DMC -providing the lowest 
possible variational energy within the same nodal surface 
of tpT- and the VMC energies clearly supports the accu- 
racy of our calculation, as well as that the N = 54 is an 
accidental closed shell, and will not be considered in the 
furthcoming analisys. 

In order to optimize the WF we use the recent method 
introduced in Reffnj. devised here in an appropriate way 
to optimize a large number of parameters during the 
AMD simulation. At each iteration time t n we com- 
pute the generalized energy gradients g n = s" 1 dE g^' IC ' 
where s is the reduced overlap matrix between the log- 
arithmic WF derivatives, appropriately regularized as in 
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Refill!. Then we use the mentioned method to minimize 
the energy in the linear space spanned by p— vectors 
9ni9n-ii ■ • • 9n~p+i corresponding to the previous p iter- 
ations (p ~ 100 for large systems). The ionic positions 
are then consistently updated according to Eq. lflU]) and, 
in order to improve the QMC stability, the corresponding 
change in the electronic variational parameters is reduced 
by a factor four. This factor clearly leads to a slowing 
down of the electron dynamic, that however turns out 
to remain very close to the Born-Oppcnheimcr one, even 
by performing only one step of WF optimization (with 
p ~ 100) each time the ionic positions are changed. On 
the other hand, as shown in the insets of Fig.©, the 
small bias due to the finite time discretization At ap- 
pears to produce only an effective change of the average 
temperature (calculated by the average ion kinetic en- 
ergy) and corresponding consistent change of the internal 
energy (see left inset). 



-0.48 



-0.49 



N=128 



N=16 




agram of Hydrogen is still under debate especially for a 
possible stable low temperature ~ A00K high pressure 
(~ 300GPa) liquid phase. We show in Fig.© the evo- 
lution of the internal energy and corresponding temper- 
ature as a function of time with the proposed AMD with 
QMC forces, starting from the bec solid at r s = 1.31. 
This solid is clearly unstable because, even at ~ 1300A" 
the internal energy decreases by a huge quantity (about 
0.02iJ per proton). Although we have not studied the 
possible stability of all other solid phases yet, the inter- 
nal energy obtained with VMC at the lowest tempera- 
ture appears slightly lower than the ground state energy 
of all plausible solid phases even when estimated by the 
DMC method [14|. This is rather remarkable because the 
application of DMC to our variational states certainly 
decreases their energies, and also that our ground state 
energy is certainly below the computed internal energy 
at finite temperature. 

In conclusion we have shown that it is possible to make 
a realistic and accurate simulation of many atoms consis- 
tently with QMC forces, for a time (~ lps) comparable 
with the present ab-initio methods based on DFT. The 
first important outcome of this calculation is that the 
bec solid structure appears clearly unstable even at low 
temperatures (N — 16 and N = 128 are consistent, sug- 
gesting small size effects in this phase), where a molec- 
ular liquid with explicit pairing correlations absent in 
the metallic solid phase, has much lower internal energy. 
This liquid phase represents either a RVB Mott insulator 
or a non conventional superconductor, stabilized only by 
the strong electron correlation. It is clear that more sys- 
tematic and accurate studies are necessary to clarify this 
novel scenario, especially on the experimental side. 
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F. Becca, M. Casula, M. Fabrizio for useful discussions, 
and the excellent stability of SP5 in CINECA. 
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FIG. 2: (color online). Evolution of the internal energy and 
temperature vs the AMD with QMC forces. At = Ao = 
1.036/s, ao = 0.7fcsTa.u.. Bottom: points represent instan- 
taneous temperatures estimated by the average kinetic en- 
ergy, lines represent the target temperatures. They should 
coincide on average for At — > 0. See right inset for a target 
temperature of 1580-K". The left inset shows the correspond- 
ing energy. The average energy, pressure and temperature in 
the last 0.5ps are -0.51319 ± 0.00003ff (-0.5127 ± 0.0001H) 
364K ± 5K (364 ± 10K ) and 335 ± 2GPa (394 ± 5GPa) for 
TV = 128 (N = 16), respectively. At each iteration the statis- 
tical noise on the total energy is ~ 5000AT. 
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